library(readr)
library(dplyr)
library(leaflet)
library(sf)
library(tidyr)
library(purrr)
library(tibble)
library(RColorBrewer)
library(gridExtra)
library(geosphere)
library(lubridate)
library(patchwork)La idea de este punto es obtener para cada zona de contaminación, las medias de intensidad, ocupacion, velocidad… de cada uno de los tres tipos de vias que se encuentran en un radio dado
trafico_final <- read_csv('/Users/pablogandia/Desktop/trafico_cluster.csv')
contaminantes <- read_csv("/Users/pablogandia/Desktop/Analisis_aire/data/raw/geoespacial/estaciones_valencia.csv") %>%
select(Estación, Latitud, Longitud, Altitud)
zonas <- c("València - Av. França",
"València - Molí del Sol", "València - Pista de Silla",
"València - Politècnic", "València Port llit antic Túria",
"València Port Moll Trans. Ponent")
contaminantes <- contaminantes[contaminantes$Estación %in% zonas, ]
datos_vias <- read_csv('/Users/pablogandia/Desktop/trafico_horario_clean.csv') %>%
select(IdPM, Hora, Dia,Mes,IntensidadMediana,VelocidadMedia,OcupacionMedia,OcupacionMaxima,DevIntensidad)Aquí establecemos para visualizar un radio de 1km, que puede ser interesante para analizar relaciones cercanas pero tampoco obviar las vias que tiene a una distancia media, que pueden ser muchas:
pal <- colorFactor(palette = brewer.pal(3, "Set1"),
domain = trafico_final$tipo_via)
m <- leaflet(trafico_final) %>%
addTiles() %>%
setView(lng = mean(trafico_final$Lon), lat = mean(trafico_final$Lat), zoom = 11) %>%
addCircleMarkers(
lng = ~Lon, lat = ~Lat,
radius = 6,
color = ~pal(tipo_via),
fillOpacity = 0.8,
stroke = FALSE,
popup = ~paste("Zona:", idpm, "<br>Tipo de vía:", tipo_via)
) %>%
addLegend(
position = "bottomright",
pal = pal,
values = ~tipo_via,
title = "Tipos de vía",
opacity = 1
) %>%
addCircles(
data = contaminantes,
lng = ~Longitud,
lat = ~Latitud,
radius = 1000,
color = "#ff9999",
fillOpacity = 0.2,
stroke = TRUE
) %>%
addCircleMarkers(
data = contaminantes,
lng = ~Longitud,
lat = ~Latitud,
radius = 8,
color = "black",
fillOpacity = 1,
popup = ~Estación
)
mVemos que principalmente las 4 estaciones que mas nos interesan son El Politécncio, La Avenida de Francia, la Pista de Silla y Moli del sol.
Las estaciones del puerto, aunque en Valencia, se ecuentran bastante distanciadas de las zonas de tráfica, y además están sesgadas por el tráfico marítimo, lo que puede complicar los analisis.
gdf_vias <- st_as_sf(trafico_final, coords = c("Lon", "Lat"), crs = 4326)
gdf_cont <- st_as_sf(contaminantes, coords = c("Longitud", "Latitud"), crs = 4326)
gdf_vias_m <- st_transform(gdf_vias, 3857)
gdf_cont_m <- st_transform(gdf_cont, 3857)
radio <- 1000
gdf_cont_buf <- st_buffer(gdf_cont_m, dist = radio)
estacion_a_vias <- map2(
gdf_cont_buf$geometry, gdf_cont$Estación,
~ {
vias_dentro <- gdf_vias_m[st_within(gdf_vias_m, .x, sparse = FALSE), ]
setNames(list(vias_dentro$idpm), .y)
}
) %>% reduce(c)
estacion_a_vias[["València - Molí del Sol"]]## [1] 1702 1715 1716 1717 1718 1719 3601 3612 3613 3622 3624 3625 3638 3639 3640
## [16] 3822vias_info <- trafico_final %>%
select(idpm, tipo_via)
df_estacion_vias <- tibble(
estacion = rep(names(estacion_a_vias), lengths(estacion_a_vias)),
idpm = unlist(estacion_a_vias)
)
df_estacion_vias <- df_estacion_vias %>%
left_join(vias_info, by = "idpm") %>%
mutate(tipo_via = factor(tipo_via, levels = c("alto", "bajo", "congestionado")))
contadores <- df_estacion_vias %>%
count(estacion, tipo_via) %>%
pivot_wider(
names_from = tipo_via,
values_from = n,
values_fill = 0,
names_prefix = "tipo_"
)
contadoresAquí vemos la cantidad de vias de diferente tipo que se le asignan a cada estacion poniendo 1km de radio.
trafico <- trafico_final %>%
left_join(datos_vias, by = c("idpm" = "IdPM"))
trafico <- trafico%>%
select(-nombre,-Lat,-Lon)
traficocontaminacion <- read_csv('/Users/pablogandia/Downloads/Calidad_Aire_limpio.csv')
contaminacion <- contaminacion %>%
mutate(
Año = year(FechaHora),
Mes = month(FechaHora),
Dia = day(FechaHora),
Hora = hour(FechaHora)
)
zonas_interes <- c("València - Av. França",
"València - Molí del Sol",
"València - Pista de Silla",
"València - Politècnic",
"València Port llit antic Túria",
"València Port Moll Trans. Ponent")
contaminacion<- contaminacion %>%
filter(Origen %in% zonas_interes, Año == 2023)generar_base_estacion <- function(nombre_estacion, trafico, contaminacion, estacion_a_vias) {
ids <- estacion_a_vias[[nombre_estacion]]
trafico_filtrado <- trafico %>%
filter(idpm %in% ids)
global <- trafico_filtrado %>%
group_by(Mes, Dia, Hora) %>%
summarise(
IntensidadGlobal = mean(IntensidadMediana, na.rm = TRUE),
VelocidadGlobal = mean(VelocidadMedia, na.rm = TRUE),
OcupacionGlobal = mean(OcupacionMaxima, na.rm = TRUE),
.groups = "drop"
)
por_tipo <- trafico_filtrado %>%
group_by(Mes, Dia, Hora, tipo_via) %>%
summarise(
Intensidad = mean(IntensidadMediana, na.rm = TRUE),
Velocidad = mean(VelocidadMedia, na.rm = TRUE),
Ocupacion = mean(OcupacionMaxima, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_wider(
names_from = tipo_via,
values_from = c(Intensidad, Velocidad, Ocupacion),
names_glue = "{.value}_{tipo_via}",
values_fill = 0
)
trafico_final <- left_join(global, por_tipo, by = c("Mes", "Dia", "Hora"))
cont_filtrada <- contaminacion %>%
filter(Origen == nombre_estacion) %>%
mutate(
Mes = month(FechaHora),
Dia = day(FechaHora),
Hora = hour(FechaHora)
)
base_final <- inner_join(cont_filtrada, trafico_final, by = c("Mes", "Dia", "Hora"))
return(base_final)
}
base_franca <- generar_base_estacion("València - Av. França", trafico, contaminacion, estacion_a_vias)
base_francaAsí quedarían los datos cruzados con contaminación (misma zona y misma fecha), y ahora además le vamos a añadir la dirección del viento de los datos del objetivo 2.
colSums(is.na(base_franca))## FechaHora Origen NO
## 0 0 0
## NO2 NOx O3
## 0 0 0
## PM10 PM2.5 SO2
## 0 0 0
## Año Mes Dia
## 0 0 0
## Hora IntensidadGlobal VelocidadGlobal
## 0 0 0
## OcupacionGlobal Intensidad_alto Intensidad_bajo
## 0 0 0
## Intensidad_congestionado Velocidad_alto Velocidad_bajo
## 0 0 0
## Velocidad_congestionado Ocupacion_alto Ocupacion_bajo
## 0 0 0
## Ocupacion_congestionado
## 0Algunos valores de cong y alta tienen nulos ya que se dan en zonas en las que no hay zonas de alta intensidad ni zonas de congestion frecuente. Por lo que estos nulos los imputaremos por 0.
viento <- read_csv('/Users/pablogandia/Downloads/meteorología_horaria_limpio.csv')
viento <- viento %>%
select(fecha,hora,"Dirección del viento")
viento <- viento %>%
mutate(
fecha = as.Date(fecha),
Mes = month(fecha),
Dia = day(fecha),
Hora = as.integer(hora)
)
simplificar_direccion <- function(dir) {
case_when(
dir %in% c("N", "NNE", "NNW") ~ "N",
dir %in% c("NE", "ENE") ~ "NE",
dir %in% c("E", "ESE") ~ "E",
dir %in% c("SE", "SSE") ~ "SE",
dir %in% c("S", "SSW") ~ "S",
dir %in% c("SW", "WSW") ~ "SW",
dir %in% c("W", "WNW") ~ "W",
dir == "NW" ~ "NW",
TRUE ~ NA_character_
)
}
viento <- viento %>%
mutate(direccion = simplificar_direccion(`Dirección del viento`))
viento <- viento %>%
select(Mes,Dia,Hora,direccion)
vientobase_franca <- base_franca %>%
left_join(viento, by = c("Mes", "Dia", "Hora")) base_franca <- base_franca %>%
mutate(
viento_Norte = as.integer(direccion %in% c("N", "NNE", "NNW")),
viento_Sur = as.integer(direccion %in% c("S", "SSW", "SSE")),
viento_Este = as.integer(direccion %in% c("E", "ESE", "ENE")),
viento_Oeste = as.integer(direccion %in% c("W", "WSW", "WNW")),
viento_Noreste = as.integer(direccion == "NE"),
viento_Sureste = as.integer(direccion == "SE"),
viento_Suroeste = as.integer(direccion == "SW"),
viento_Noroeste = as.integer(direccion == "NW"),
viento_Sin_dato = as.integer(is.na(direccion))
)
base_francaEl resultado final es una base de datos que estudia contaminación en una zona y fecha concreta, asi como la direccion del viento y el tráfico de sus alrededores bastante bien carcterizado y segmentado.